!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the subroutine is used to calcaulate interior point 
! papramenters by MOC.
! ref. to chapter 16 of 'Gas Dynamics' by Maurice J. Zucrow & Joe D. Hoffmann
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-06-16
! revised  :
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! include  :
!   subroutine : InterP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   x,y       point location, m
!   Vu,Vv     axial, radical velocity, m/s
!   icor      correction steps
!   i         cycle index
!   theta     streamline angle, rad
!   a         speed of sound, m/s
!   al        alpha, Mach angle, rad
!   l,Q,R,S,T equation coefficients
!
!   subscript:
!      p   plus(except Cp), +
!      m   minus, -
!      f   front cycle    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   1. angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine InterP()
   use VariableDef
   implicit none
!         1       
!           4()
!       2
!       given variable
!   real*8::x1,y1,Vu1,Vv1,x2,y2,Vu2,Vv2
!       temp variable
   real*8::Vup,Vvp,yp,thetap,Vp,ap,Mp,alp,Vum,Vvm,ym,thetam,Vm,am,Mm,alm
   real*8::lp,Qp,Rp,Sp,Tp,lm,Qm,Rm,Sm,Tm
   real*8::x4f,y4f,Vu4f,Vv4f,a4
!       calculate variable
!   real*8::x4,y4,Vu4,Vv4,V4,p4,T4,rou4,M4
!       cycle index
   integer::icor,i
  
   Vup=Vu2
   Vvp=Vv2
   yp=y2
   Vum=Vu1
   Vvm=Vv1
   ym=y1
   x4f=0.0
   y4f=0.0
   Vu4f=0.0
   Vv4f=0.0
   icor=6

!   call GasCons((T1+T2)/2,Rg,g,Cp)

   do i=1,icor
      thetap=atan(Vvp/Vup)
      Vp=sqrt(Vup**2+Vvp**2)
      call Thermo(Vp,p4,rou4,T4,ap,Mp)
      alp=asin(1/Mp)
             
      thetam=atan(Vvm/Vum)
      Vm=sqrt(Vum**2+Vvm**2)
      call Thermo(Vm,p4,rou4,T4,am,Mm)
      alm=asin(1/Mm)
      
      lp=tan(thetap+alp)
      Qp=Vup**2-ap**2
      Rp=2*Vup*Vvp-Qp*lp
!    corrected for point on the axis
      if (yp==0.0) then
         Sp=delta*Vvm*ap**2/ym
      else
         Sp=delta*Vvp*ap**2/yp
      end if
!    end corrected for point on the axis
             
      lm=tan(thetam-alm)
      Qm=Vum**2-am**2
      Rm=2*Vum*Vvm-Qm*lm
      Sm=delta*Vvm*am**2/ym
      
      x4=(lp*x2-lm*x1+y1-y2)/(lp-lm)
      y4=(lp*lm*(x2-x1)+lp*y1-lm*y2)/(lp-lm)
      
      if (y4<0.0) then   ! y4<0, exit
         flagy40=1
         return
      end if
      
      Tp=Sp*(x4-x2)+Qp*Vu2+Rp*Vv2
      Tm=Sm*(x4-x1)+Qm*Vu1+Rm*Vv1
      
      Vu4=(Tp*Rm-Tm*Rp)/(Qp*Rm-Qm*Rp)
      Vv4=(Tm*Qp-Tp*Qm)/(Qp*Rm-Qm*Rp)
      
      if ((abs(x4-x4f) <= el) .or. (abs(y4-y4f) <= el)) then
         if ((abs(Vu4-Vu4f) <= ev) .and. (abs(Vv4-Vv4f) <= ev)) then
            exit
         end if
      end if
             
      Vup=(Vu2+Vu4)/2
      Vvp=(Vv2+Vv4)/2
      yp=(y2+y4)/2
      Vum=(Vu1+Vu4)/2
      Vvm=(Vv1+Vv4)/2
      ym=(y1+y4)/2

      x4f=x4
      y4f=y4
      Vu4f=Vu4
      Vv4f=Vv4 
   end do
   V4=sqrt(Vu4**2+Vv4**2)
   call Thermo(V4,p4,rou4,T4,a4,M4)
        
end subroutine InterP
